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Abstract 

We used a molecular dynamics simulation with the modified Brenner reactive empirical bond order potential to 
investigate the erosion of a graphite surface due to the incidence of hydrogen, deuterium, and tritium atoms. Incident 
particles cause pressure on the graphite surface, and the chemical bond between graphene layers then generates 
heat to erode the graphite surface. We evaluated the speed of surface destruction by calculating the pseudo-radial 
distribution function. The speed of surface destruction due to incident hydrogen isotopes was higher than that due to 
hydrogen atoms. The surface destruction increased exponentially and its decay time constant was a power function 
of the incident energy. We measured the erosion yield, which indicated a steady state for the graphite erosion. The 
erosion yield flux in the steady state increased linearly with the incident energy. The erosion yield flux was almost 
independent of the type of incident particle, and the erosion yield start time was smaller for hydrogen isotopes than 
for hydrogen atoms. 
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1. Introduction 

In the context of research into nuclear fusion, 
we studied the plasma surface interaction (PSI) 
problem [1,2,3,4,5]. A portion of the plasma con- 
fined in an experimental device falls onto a diverter 
wall, which is a shield made of graphite or carbon 
fiber composite tiles. The incident hydrogen plasma 
erodes these carbon tiles in a process called chem- 
ical sputtering. The erosion produces hydrocarbon 
molecules, such as CH X and C2rl x , which affect the 
plasma confinement. 
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To solve the PSI problem, the mechanism of the 
graphite erosion has been researched using molecu- 
lar dynamics simulation (MD) [6,7,8,9]. Previously, 
we investigated the PSI of graphite surfaces using 
the modified Brenner reactive empirical bond or- 
der (REBO) potential [10]. That MD simulation 
showed that if incident energy was 5 eV, almost 
of all incident hydrogen atoms were absorbed by 
the graphite surface, while if the incident energy 
was 15 eV, most incident hydrogen atoms were re- 
flected. This absorption and reflection can be ex- 
plained by the chemical reaction between a single hy- 
drogen atom and a single graphene [11,12]. However, 
the number of absorbed hydrogen atoms seems to 
be independent of this graphite erosion because al- 
though the hydrogen atoms are absorbed by the first 
graphene on the surface only, multiple graphencs arc 
destroyed simultaneously. Almost all of the absorbed 
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Fig. 1. The erosion of the graphite surface for deuterium incident at 5 eV. 



hydrogen atoms are located on one side of the first 
graphene. The chemical bond between the first and 
second graphcnes triggers the graphite erosion. We 
assume that the momentum of the incident hydro- 
gen atoms, which is considered to be pressure from 
the point of view of macroscopic thermodynamics, 
presses the first graphene against the second. 

Nuclear fusion research must not only consider 
the hydrogen atom, but also the hydrogen isotopes 
deuterium and tritium. It is important to under- 
stand the difference in the PSI due to these isotopes. 
A study of the chemical reaction between a single 
graphene and a single hydrogen isotope atom showed 
that the absorption and reflection rates as functions 
of the incident energy differed from that of the hy- 
drogen atom [13]. These differences will be also be 
reflected in the graphite erosion. Moreover, the in- 
cident momentum of the hydrogen isotopes differs 
than that of the hydrogen atoms even if the incident 
energies are the same. Therefore, the pressure of hy- 
drogen isotopes on the graphite surface is greater 
than that of hydrogen atoms. If the chemical bond- 
ing between the first and second graphcnes, caused 
by pressure from the incident momentum, triggers 



the graphite erosion, we would expect that the hy- 
drogen isotopes would accelerate this erosion more 
than hydrogen atoms would. 

We used MD simulation to investigate graphite 
erosion due to the incidence of hydrogen, deuterium, 
and tritium atoms. We describe the simulation 
model and method in §2. In §3, we present and dis- 
cuss the simulation results. This paper concludes 
with a §4. 

2. Simulation Method 

As the graphite, we located eight graphcnes [14] 
with an "ABAB" lattice structure parallel to x-y 
plane. Each graphene consisted of 160 carbon atoms 
measuring 2.13 nm x 1.97 nm. The size of the sim- 
ulation box in the x— and y-directions equaled that 
of graphene with the periodic boundary condition. 
The inter-layer distance of the graphite was initially 
3.35 A. The initial distribution of the momentum of 
the carbon atoms followed the Maxwell-Boltzmann 
distribution at 300 K. During the simulation, only 
two carbon atoms were fixed to support the base of 
the graphite. One was the center atom of the 7-th 
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Fig. 2. The erosion of the graphite surface for deuterium incident at 15 eV. 



graphcne from the surface, and the other was located 
at the boundary of the 8-th graphene. The graphite 
surface was oriented to face the positive ^-direction. 

For the simulation, 500 or more hydrogen or iso- 
tope atoms were injected at regular time intervals 
of 0.1 ps. The x— and y-coordinates of the injec- 
tion point were set at random. The z-coordinate of 
the injection point was 60 A. The initial momentum 
vector (0, 0, po) was parallel to the z-axis, and was 
defined by 



i,j>i 



(2) 



where is the distance between the i— th and j-th 
atoms. The functions Vj^j and represent repul- 
sion and attraction, respectively. The function bij 
generates multi-body force. To conserve the accu- 
racy of the calculation, the time step was 5 x 10 -18 s. 
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3. Results and Discussion 



where E\ is the incident energy, and m is the mass of 
the incident particle, which is 1. 2. or 3 u for the case 
of hydrogen, deuterium, or tritium, respectively. 

We performed our MD simulation under NVE 
conditions, where the number of atoms, volume, and 
total energy are conserved, except for the addition 
of incident atoms and removal of outgoing atoms. 
The simulation time was developed using second 
order symplcctic integration [15]. The chemical in- 
teraction was represented by the modified Brenner 
REBO potential [12,16]: 



We performed simulations for incident energies 
in the range 0.5-30 eV. Figure 1 shows a snap- 
shot of the graphite erosion due to deuterium in- 
cidence at 5 eV. This figure clearly illustrates the 
erosion process. In the initial short time period (5 
ps), the deuterium atoms are generally absorbed. 
This is consistent with previous research findings 
concerning the chemical reaction between a single 
deuterium atom and a single graphene. Although 
the first graphene absorbs many deuterium atoms, 
no erosion yet occurs. When the first and second 



graphenes are linked by a covalent bond (12 ps), 
the erosion process starts, using the binding en- 
ergy of the chemical bond between the graphenes. 
The graphite erosion then advances to the lower 
graphenes connected by covalent bonds (28 ps). 
For deuterium incident at 15 eV (see Fig. 2), al- 
though the first graphene reflects almost all of the 
deuterium atoms (2 ps), and the chemical bonding 
between the graphenes results in graphite erosion (5 
ps), as in the case of 5 eV. This behavior is the same 
for both hydrogen atoms and hydrogen isotopes. 
Therefore, we conclude that the chemical bonding 
between the graphenes is the dominant mechanism 
in the erosion process, rather than the absorption 
and reflection at the graphene surface. 



a) Deuterium at 5 eV 
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b) Tritium at 5 eV 
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To estimate the surface destruction, we must 
first calculate the radial distribution function. How- 
ever, we cannot define the volume in this simula- 
tion model because there is no boundary in the 
z-direction. Therefore, we use a pseudo-radial dis- 
tribution function g(r,t) defined by 



9(r,t) 



1 dn(r, t) 
47rr 2 dr 



(3) 



where r is the distance between two particles, t is 
time, and n(r, t) is the number of carbon atoms lo- 
cated at a distance of less than r at time t. Although 
the pseudo-radial distribution function g(r, t) differs 
from usual radial distribution function in that it is 
not normalized by the number density, it is sufficient 
to estimate the lattice structure because the num- 
ber density is constant. Figure 3 shows the change of 
the pseudo-radial distribution function g(r, t) with 
time. It is the same for hydrogen, deuterium, and 
tritium. Initially, the peaks of the graphite lattice 
structure are sharp, but then they broaden as time 
goes on, indicating the progress of the surface de- 
struction. To measure the speed of surface destruc- 
tion, we plotted the maximum values of the pseudo- 
radial distribution function g m ax(fjt) as a function 
of time (see Fig. 4). These values often indicate C- 
C bonds of length r = 1.42 Aand correspond to the 
number of sp 2 bonds. From this figure, if we neglect 
the irregular regions, which are for t > 35 ps at 20 
eV and for t > 25 ps at 30 eV, it is obvious that the 
number of sp 2 bonds decreases exponentially with 
time. The speed of the surface destruction due to 
incident hydrogen isotopes is greater than that due 
to hydrogen atom. Figure 5 shows that the decay 
time constant r of the surface destruction is a power 
function of the incident energy E\. We found that 
the maximum value of the pseudo-radial distribu- 
tion function g max (r, t) and the decay time constant 
r can be represented by 



(r, t) cx exp 



t 



t = CEj a , 



(4) 
(5) 



Fig. 3. The pseudo— radial distribution function g(r, t) as a 
function of the distance between carbon atoms r and time t. 



We investigated the surface destruction (amor- 
phization of the graphite surface) and the erosion 
yield in our study of graphite erosion. 



where the scaling exponent a for the hydrogen, deu- 
terium, and tritium atoms is —0.792, —0.778, and 
—0.767, respectively. The constant value C for the 
hydrogen, deuterium, and tritium atoms is 1.37 x 
lO" 10 , 1.12 x 10~ 10 , and 1.03 x lO" 10 , respectively. 
This indicates that the scaling exponent a is almost 
independent of the type of incident particle, while 
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a) 2 eV 
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Fig. 4. The maximum value of the pseudo-radial distribution function g m ax(r, t) as a function of the time t. 



the constant value C does depend on the incident 
particle type. 
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Fig. 5. The decay time constant r as a function of the inci- 
dent energy. The solid, long dashed, and short dashed lines 
indicate values determined by the simulation for incident 
hydrogen, deuterium, and tritium, respectively. 

Next, we observe the erosion yield. In our simu- 
lation, almost all of the molecules produced had a 
chain structure. That is to say, carbon atoms were 
bounded by sp 1 bonds, and hydrogen atoms and iso- 
topes often terminated the end of a carbon chain 
(e.g., H-C-C-C-C-H) . Since a chemical reaction 
occurs between these molecules and changes their 
molecular structure, we did not specify the types 
of molecule produced. For the erosion yields Y, we 
counted the number of carbon atoms that moved to 
the region z > 24 A, where the first graphene is ini- 
tially located on z = 11.7 A. Fig. 6 shows the ero- 
sion yield Y as a function of time t. First, the yield 
molecules are created after an initial delay. This im- 
plies that chemical sputtering requires an increase in 
surface temperature. Second, the erosion yield Y in- 
creases linearly to a level of 740 carbon atoms. This 
linear increase of Y is thought to be the graphite ero- 
sion steady state. After reaching 740 carbon atoms, 
the erosion yield Y increases rapidly to the point of 
saturation. It is clear from our MD simulations that 
the base graphene is destroyed when Y > 740 car- 
bon atoms. The eroded graphene appears to trans- 
fer heat to lower graphenes during the steady state, 
and the rapid increase and saturation (Y > 740) are 
due to the absence of additional lower graphenes to 
absorb the heat. We fit the erosion yields Y to the 
following linear function yi(t): 

yi (t) = ^S(t-t ), (6) 



where the fitted region is 5 < Y < 740. The coeffi- 
cient S is the square measure of 2.13 nm x 1.97 nm. 
The fitting parameters <j) and to correspond to the 
erosion yield flux and the erosion yield start time, 
respectively. Fig. 7 shows the incident energy and 
isotope dependence of the erosion yield flux. The 
erosion yield flux increases as a linear function of the 
incident energy, but it is almost independent of the 
type of isotope. It is clear that in the steady state, the 
erosion yield depends on the incident energy rather 
than the incident momentum, which depends on the 
type of incident particle. Moreover, the linear in- 
crease of the erosion yield flux indicates that the to- 
tal yield on a real graphite surface also increases lin- 
early with the incident energy. In other experiments 
[17,18] , the total yield increased linearly with the in- 
cident energy up to 50 eV. Fig. 8 shows the incident 
energy and isotope dependence of the erosion yield 
start time. This indicates that incident hydrogen iso- 
topes start the erosion more quickly than hydrogen 
atom does. In addition, the erosion yield start time 
decreases as the incident energy increases, except 
with tritium at 2 eV. Our explanation of this is as 
follows. The erosion yield requires an increase in the 
graphite surface temperature, and it is the chemi- 
cal bonding between the first and second graphenes 
that generates the heat to cause this surface tem- 
perature increase. Because the graphenes are linked 
due to a pressure on surface, which is derived from 
incident momentum, the erosion yield start time de- 
pends on the type of incident particle. 

Next, we address two other important problems. 

The first is the problem of graphite interlayer in- 
teraction. In this work, although the interaction of a 
covalcnt bond is represented by the modified Bren- 
ner REBO potential model, the graphite interlayer 
intermolecular force is not included. While the inter- 
action energy between graphite layers is thought to 
be about 0.01 times that of a covalent bond, its true 
value is not yet known [19]. In addition, when us- 
ing simple two-body interaction functions, such as 
the van der Waals interaction in the graphite inter- 
layer interaction model, the real graphite "ABAB" 
structure cannot be represented. Therefore, we do 
not have an intermolecular interaction model that is 
suitable for graphite. However, the modified Bren- 
ner REBO potential generates a strong repulsive 
force when the interlayer distance is less than 2 A. 
Although the interlayer distance cannot be kept at 
3.35 A, the repulsive force maintains the layer struc- 
ture of the graphite as long as the graphite is not 
subjected to a strong external force or high pres- 
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Fig. 6. The erosion yields Y of carbon atoms as a function of time. 



sure, and its temperature is maintained below the 
melting point. In his paper on the chemical sput- 
tering of polymer, Yamashiro described the impor- 
tance of the intcrmolccular attractive force in cre- 
ating a cluster structure and bringing about ther- 
mal deposition [20]. However, the intermolecular re- 
pulsive force, which maintains the interlayer dis- 
tance and "ABAB" structure, plays a more impor- 
tant role than the intermolecular attractive force 
in the graphite structure. In addition, the hydro- 
carbon molecules that are produced tend to prefer 
chain structures to cluster structures. Therefore, we 
determined that it was necessary to investigate the 
importance of the interlayer intermolecular force in 
the specific case of graphite erosion. We have been 
creating a new intermolecular potential model for 
graphite interlayer interaction. This work is essen- 
tial for comparison with our future work, which will 
include the graphite interlayer intermolecular inter- 
action. 

The second problem relates to temperature con- 
trol. In this work, we can achieve steady state 
graphite erosion without the need to cool the 
graphite, since heat is conducted to the lower 
graphencs. In other words, the lower graphenes act 
as a heat reservoir. MD simulations generally use a 
thermostat to provide rapid cooling because of the 
picosecond-scale simulation time. Consequently, 
the dynamics of atoms in the simulation differ from 
those in reality; for example, they halt the reac- 
tion due to a trap on the local potential minimum. 
Moreover, because the rate of heat conduction de- 
pends on the interaction potential model among 
atoms, which represents a realistic interaction, it is 
fairly realistic but slow, judging from the picosec- 
ond simulation time scale, and cannot be changed 
too quickly. Therefore, even if the thermostat were 
not set to the position where a chemical reaction 
occurs simply to preserve the realistic dynamics of 
atoms, sufficient cooling would still not be possible 
because of the slow rate of heat conduction. There- 
fore, we believe that the steady state of the graphite 
erosion that we have identified here is a realistic 
phenomenon, even if the steady state 1 is maintained 
for a short time only, and the incident flux is very 
high. A longer-term steady state would likely be 
possible using graphite with more graphene layers. 
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Fig. 7. The erosion yield flux <j> of carbon atoms as a func- 
tion of the incident energy. The solid, long dashed, and short 
dashed lines indicate the flux due to hydrogen, incident deu- 
terium, and tritium, respectively. 
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4. Summary 

We have studied the graphite erosion process due 
to the incidence of hydrogen, deuterium, and tri- 
tium atoms using MD simulation with the modi- 
fied Brenner REBO potential for incident energy Ej 
in the range 0.5-30 eV. For E\ = 5 eV, both hy- 
drogen atoms and hydrogen isotopes were generally 
absorbed by the graphene surface. At E\ — 15 eV, 
however, almost all of the incident particles were re- 
flected. When the incident particles press the first 
graphene against the second, erosion of the graphite 
surface is caused by the binding energy of the chemi- 
cal bond between the graphene layers. We evaluated 
the speed of the surface destruction by calculating 
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the pseudo-radial distribution function g{r,t) 1 the 
peak values of which decrease with time. The max- 
imum value of the pseudo-radial distribution func- 
tion j max (r, t) . which generally indicates the number 
of sp 2 C-C bonds, decreased exponentially. The sur- 
face destruction due to the hydrogen isotopes pro- 
gressed faster than the destruction caused by hy- 
drogen atoms. Moreover, we found a power law be- 
tween the decay time constant and incident energy, 
and the scaling exponent seems to be independent 
of the type of incident particle. Measuring the ero- 
sion yields indicated a steady state of the graphite 
erosion. The erosion yield flux in the steady state 
increased linearly with the incident energy, which 
agreed with experimental results. The erosion yield 
flux is almost independent of the type of incident 
isotope. Observation of the erosion yield start time 
showed that the erosion due to hydrogen isotopes 
started earlier than the erosion due to the hydrogen 
atoms. 
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